KELPER2

Dataset EDA

This document is for preliminary exploration of the data from the Subtidal Ecology lab at Pontificia Universidad. See data-collection.html for details on the collection.

1 Allometry

Start with the simplest. Holdfast diameter will be the primary metric for tracking plant size and will be used to estimate everything else, including biomass and demographic rates.

morph.df <- list(
  IFOP_morph=read_csv(glue("{data.dir}IFOP_Kelp_Morphology.csv"),
                      na=c("NaN", "")) %>% 
    select(c(1:9,11:13,15:19)) %>%
    filter(Group=="Macroalga" &
             Species=="Lessonia trabeculata" &
             Abundance > 0) %>% 
    mutate(Date=ymd(paste(Year, Month, Day, sep="-")),
           Site=factor(Site),
           Management=factor(Management),
           Frondosity=factor(Frondosity),
           source="IFOP_morph") %>%
    rename(LengthTot=TL),
  IFOP_allom=read_csv(glue("{data.dir}IFOP_allometry.csv")) %>% 
    select(-16) %>%
    mutate(Date=ymd(paste(Year, Month, Day, sep="-")),
           Site=factor(Site),
           Management=factor(Management),
           Weight=1e3*Weight,
           source="IFOP_allom") %>%
    rename(LengthTot=TL,
           WeightTot=Weight),
  NERC_bl=read_csv(glue("{data.dir}Quads_baseline.csv")) %>% 
    select(-c(1,3:4,10:12)) %>%
    filter(Species=="Lessonia trabeculata") %>%
    mutate(Site=factor(Site),
           Zone=factor(Zone),
           Management=factor(Management),
           Upwelling=factor(Upwelling),
           Frondosity=factor(Frondosity),
           source="NERC_bl") %>%
    rename(Holdfast=Holdfast.Diameter,
           LengthTot=Total.Length),
  NERC_ph=read_csv(glue("{data.dir}Patches_Morphology.csv")) %>%
    mutate(Site=factor(Site),
           Zone=factor(Zone),
           Management=factor(Management),
           Upwelling=factor(Upwelling),
           source="NERC_ph") %>%
    rename(Holdfast=Holdfast.Diameter,
           LengthTot=Total.Length),
  Cata_ph=read_csv(glue("{data.dir}Morfo_DW_C.csv")) %>%
    select(-c(1,11,20:30)) %>%
    mutate(Site=factor(Site),
         Zone=factor(Zone),
         Management=factor(Management)) %>%
    rename(Holdfast=Diam_Holdfast,
           Height=High_Holdfast,
           LengthTot=LT,
           WeightTot=TW_sum,
           Weight=Weight_Str,
           Length=LT_stipes,
           Dichot=Dichot_stipes,
           Diam=Diam_stipes) %>%
    pivot_wider(names_from="Structure", values_from=c(8:17)) %>%
    select(!where(~all(is.na(.x)))) %>%
    select(-starts_with("Structure_")) %>%
    mutate(source="Cata_ph") %>%
    rename(Stipes=Stipes_Stipes)
) %>%
  do.call('bind_rows', .)

glimpse(morph.df)

## Rows: 4,563
## Columns: 49
## $ Serie              <dbl> 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 1…
## $ Date               <date> 2018-05-16, 2018-05-16, 2018-05-16, 2018-05-16, 20…
## $ Day                <dbl> 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16,…
## $ Month              <dbl> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, …
## $ Year               <dbl> 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 201…
## $ Site               <fct> Quintay A, Quintay A, Quintay A, Quintay A, Quintay…
## $ Management         <fct> TURF, TURF, TURF, TURF, TURF, TURF, TURF, TURF, TUR…
## $ Transect           <dbl> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, …
## $ Quadrats           <dbl> 1, 2, 2, 3, 4, 5, 5, 6, 6, 6, 6, 6, 6, 6, 7, 7, 8, …
## $ Depth              <dbl> 20.0, 18.8, 18.8, 19.0, 18.7, 18.8, 18.8, NA, NA, N…
## $ Group              <chr> "Macroalga", "Macroalga", "Macroalga", "Macroalga",…
## $ Species            <chr> "Lessonia trabeculata", "Lessonia trabeculata", "Le…
## $ Abundance          <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ Holdfast           <dbl> 33, 16, 25, 23, 23, 37, 38, 25, 10, 7, 7, 22, 12, 2…
## $ Stipes             <dbl> 16, 4, 3, 6, 7, 7, 7, 7, 5, 3, 14, 12, 3, 15, 10, 1…
## $ LengthTot          <dbl> 128, 179, 133, 160, 118, 171, 213, 60, 50, 100, 100…
## $ Frondosity         <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ source             <chr> "IFOP_morph", "IFOP_morph", "IFOP_morph", "IFOP_mor…
## $ Extraction         <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ WeightTot          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Holdfast_perimeter <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Reproductive       <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Zone               <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Upwelling          <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Station            <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Quad               <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ ...1               <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Survey.Number      <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ MAE                <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Latitude           <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Longitude          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Surveyor.Name      <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Patch              <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Replicate          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ ID                 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Weight_Holdfast    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Weight_Stipes      <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Weight_Blades      <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ FW_Holdfast        <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ FW_Stipes          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ FW_Blades          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ DW_Holdfast        <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ DW_Stipes          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ DW_Blades          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Height_Holdfast    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ First_bif_Stipes   <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Diam_Stipes        <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Length_Stipes      <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Dichot_Stipes      <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…

summary(morph.df)

##      Serie             Date                 Day            Month      
##  Min.   :   1.0   Min.   :2018-04-20   Min.   : 2.00   Min.   : 1.00  
##  1st Qu.: 993.8   1st Qu.:2019-11-10   1st Qu.: 7.00   1st Qu.: 4.00  
##  Median :2618.5   Median :2020-10-04   Median :16.00   Median : 8.00  
##  Mean   :2257.7   Mean   :2020-06-12   Mean   :14.91   Mean   : 7.22  
##  3rd Qu.:3407.2   3rd Qu.:2021-04-19   3rd Qu.:24.00   3rd Qu.:11.00  
##  Max.   :4240.0   Max.   :2022-05-17   Max.   :30.00   Max.   :12.00  
##  NA's   :2131     NA's   :87           NA's   :2131    NA's   :2131   
##       Year                             Site        Management  
##  Min.   :2018   Isla Grande de Atacama   : 503   OA     :1037  
##  1st Qu.:2018   La Ballena               : 382   Reserve: 154  
##  Median :2020   Chanaral de Aceituno OA  : 370   TURF   :3027  
##  Mean   :2020   Puerto Viejo             : 311   AMERB  : 345  
##  3rd Qu.:2021   Cobquecura               : 287                 
##  Max.   :2022   Chanaral de Aceituno TURF: 280                 
##  NA's   :2131   (Other)                  :2430                 
##     Transect         Quadrats          Depth          Group          
##  Min.   : 1.000   Min.   : 1.000   Min.   : 3.30   Length:4563       
##  1st Qu.: 2.000   1st Qu.: 3.000   1st Qu.: 9.10   Class :character  
##  Median : 3.000   Median : 4.000   Median :11.00   Mode  :character  
##  Mean   : 4.167   Mean   : 4.537   Mean   :11.21                     
##  3rd Qu.: 6.000   3rd Qu.: 6.000   3rd Qu.:12.80                     
##  Max.   :10.000   Max.   :10.000   Max.   :24.60                     
##  NA's   :1723     NA's   :2451     NA's   :1806                      
##    Species            Abundance         Holdfast         Stipes      
##  Length:4563        Min.   :0.0000   Min.   : 0.00   Min.   : 0.000  
##  Class :character   1st Qu.:1.0000   1st Qu.: 7.00   1st Qu.: 2.000  
##  Mode  :character   Median :1.0000   Median :19.00   Median : 3.000  
##                     Mean   :0.9996   Mean   :19.28   Mean   : 5.091  
##                     3rd Qu.:1.0000   3rd Qu.:29.00   3rd Qu.: 6.667  
##                     Max.   :1.0000   Max.   :95.00   Max.   :79.000  
##                     NA's   :1723     NA's   :1       NA's   :86      
##    LengthTot     Frondosity     source           Extraction       
##  Min.   :  0.0   1   :  56   Length:4563        Length:4563       
##  1st Qu.: 50.0   2   : 358   Class :character   Class :character  
##  Median :131.0   3   :1621   Mode  :character   Mode  :character  
##  Mean   :126.6   0   :   1                                        
##  3rd Qu.:193.0   NA's:2527                                        
##  Max.   :450.0                                                    
##  NA's   :18                                                       
##    WeightTot       Holdfast_perimeter  Reproductive              Zone     
##  Min.   :    1.0   Min.   :  0.80     Min.   :0.000   Central      : 606  
##  1st Qu.:  332.5   1st Qu.:  9.00     1st Qu.:0.000   North        : 814  
##  Median : 1380.0   Median : 18.00     Median :0.000   North-Central: 650  
##  Mean   : 2637.0   Mean   : 20.23     Mean   :0.424   Center       :  28  
##  3rd Qu.: 3936.5   3rd Qu.: 27.50     3rd Qu.:1.000   North-Center :  33  
##  Max.   :20000.0   Max.   :161.00     Max.   :1.000   NA's         :2432  
##  NA's   :4156      NA's   :4376       NA's   :4280                        
##         Upwelling       Station            Quad             ...1       
##  nonupwelling: 503   Min.   :  0.00   Min.   : 1.000   Min.   :   1.0  
##  upwelling   : 225   1st Qu.: 20.00   1st Qu.: 3.000   1st Qu.: 329.8  
##  N           : 891   Median : 40.00   Median : 5.000   Median : 658.5  
##  Y           : 425   Mean   : 41.05   Mean   : 4.919   Mean   : 658.5  
##  NA's        :2519   3rd Qu.: 60.00   3rd Qu.: 7.000   3rd Qu.: 987.2  
##                      Max.   :100.00   Max.   :10.000   Max.   :1316.0  
##                      NA's   :3835     NA's   :3835     NA's   :3247    
##  Survey.Number        MAE           Latitude        Longitude     
##  Min.   :0.000   Min.   : 0.00   Min.   :-32.31   Min.   :-71.53  
##  1st Qu.:1.000   1st Qu.:12.00   1st Qu.:-32.21   1st Qu.:-71.50  
##  Median :3.000   Median :24.00   Median :-29.05   Median :-71.48  
##  Mean   :2.337   Mean   :19.22   Mean   :-29.41   Mean   :-71.30  
##  3rd Qu.:3.000   3rd Qu.:24.00   3rd Qu.:-27.34   3rd Qu.:-70.98  
##  Max.   :4.000   Max.   :30.00   Max.   :-27.25   Max.   :-70.96  
##  NA's   :3247    NA's   :3247    NA's   :3257     NA's   :3257    
##  Surveyor.Name          Patch         Replicate            ID        
##  Length:4563        Min.   :1.000   Min.   : 1.000   Min.   : 1.000  
##  Class :character   1st Qu.:1.000   1st Qu.: 4.000   1st Qu.: 4.000  
##  Mode  :character   Median :2.000   Median : 7.000   Median : 8.000  
##                     Mean   :2.438   Mean   : 8.362   Mean   : 8.138  
##                     3rd Qu.:3.000   3rd Qu.:12.000   3rd Qu.:11.000  
##                     Max.   :4.000   Max.   :26.000   Max.   :22.000  
##                     NA's   :3247    NA's   :3252     NA's   :4476    
##  Weight_Holdfast  Weight_Stipes    Weight_Blades   FW_Holdfast    
##  Min.   :   1.0   Min.   :   0.0   Min.   :  14   Min.   : 0.179  
##  1st Qu.:  46.0   1st Qu.:  53.0   1st Qu.: 187   1st Qu.: 1.103  
##  Median : 178.0   Median : 377.0   Median : 471   Median : 1.992  
##  Mean   : 390.2   Mean   : 831.4   Mean   :1029   Mean   : 2.595  
##  3rd Qu.: 575.5   3rd Qu.:1152.0   3rd Qu.:1088   3rd Qu.: 3.245  
##  Max.   :3200.0   Max.   :8615.0   Max.   :6885   Max.   :13.533  
##  NA's   :4476     NA's   :4476     NA's   :4476   NA's   :4476    
##    FW_Stipes       FW_Blades      DW_Holdfast      DW_Stipes    
##  Min.   :0.191   Min.   :0.023   Min.   :0.040   Min.   :0.045  
##  1st Qu.:1.050   1st Qu.:0.526   1st Qu.:0.289   1st Qu.:0.224  
##  Median :1.891   Median :0.666   Median :0.492   Median :0.485  
##  Mean   :2.311   Mean   :0.662   Mean   :0.784   Mean   :0.673  
##  3rd Qu.:3.055   3rd Qu.:0.821   3rd Qu.:0.968   3rd Qu.:0.792  
##  Max.   :7.161   Max.   :1.529   Max.   :4.041   Max.   :3.332  
##  NA's   :4477    NA's   :4476    NA's   :4476    NA's   :4477   
##    DW_Blades     Height_Holdfast  First_bif_Stipes  Diam_Stipes    
##  Min.   :0.030   Min.   :  1.00   Min.   : 0.50    Min.   :  1.00  
##  1st Qu.:0.107   1st Qu.:  5.00   1st Qu.: 2.65    1st Qu.: 12.17  
##  Median :0.132   Median : 11.00   Median : 8.10    Median : 26.60  
##  Mean   :0.140   Mean   : 12.74   Mean   :15.93    Mean   : 27.97  
##  3rd Qu.:0.177   3rd Qu.: 17.00   3rd Qu.:26.17    3rd Qu.: 37.20  
##  Max.   :0.446   Max.   :101.00   Max.   :79.00    Max.   :296.00  
##  NA's   :4476    NA's   :4476     NA's   :4476     NA's   :4478    
##  Length_Stipes     Dichot_Stipes    
##  Min.   :  0.933   Min.   :  0.933  
##  1st Qu.: 19.175   1st Qu.: 19.175  
##  Median : 62.500   Median : 62.500  
##  Mean   : 63.080   Mean   : 63.080  
##  3rd Qu.: 97.500   3rd Qu.: 97.500  
##  Max.   :187.333   Max.   :187.333  
##  NA's   :4479      NA's   :4479

ggplot(morph.df, aes(Holdfast)) + geom_histogram() + 
  facet_wrap(~source, scales="free_y")

ggplot(morph.df, aes(Holdfast, colour=source)) + geom_density()

ggplot(morph.df, aes(LengthTot)) + geom_histogram() + 
  facet_wrap(~source, scales="free_y")

ggplot(morph.df, aes(LengthTot, colour=source)) + geom_density()

ggplot(morph.df, aes(WeightTot)) + geom_histogram() + 
  facet_wrap(~source, scales="free_y")

ggplot(morph.df, aes(WeightTot, colour=source)) + geom_density()

ggplot(morph.df, aes(Frondosity)) + geom_bar() + 
  facet_wrap(~source, scales="free_y")

ggplot(morph.df, aes(factor(Reproductive))) + geom_bar() + 
  facet_wrap(~source, scales="free_y")

ggplot(morph.df, aes(Holdfast, Stipes, colour=source)) + 
  geom_point(alpha=0.5, shape=1) + 
  stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site)) +
  stat_smooth(method="lm", se=F, size=0.75, aes(group=source)) +
  scale_x_continuous("Holdfast diameter", trans="log1p") +
  scale_y_continuous(trans="log1p")

ggplot(morph.df, aes(Holdfast, LengthTot, colour=source)) + 
  geom_point(alpha=0.5, shape=1) + 
  stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site),
              formula=y~x+I(x^2)) +
  stat_smooth(method="lm", se=F, size=0.75, aes(group=source),
              formula=y~x+I(x^2)) +
  scale_x_continuous("Holdfast diameter", trans="log1p") +
  scale_y_continuous(trans="log1p")

ggplot(morph.df, aes(Holdfast, WeightTot, colour=source)) + 
  geom_point(alpha=0.5, shape=1) + 
  stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site),
              formula=y~x+I(x^2)) +
  stat_smooth(method="lm", se=F, size=0.75, aes(group=source),
              formula=y~x+I(x^2)) +
  scale_x_continuous("Holdfast diameter", trans="log1p") +
  scale_y_continuous(trans="log1p")

ggplot(morph.df, aes(LengthTot, WeightTot, colour=source)) + 
  geom_point(alpha=0.5, shape=1) + 
  stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site)) +
  stat_smooth(method="lm", se=F, size=0.75, aes(group=source)) +
  scale_x_continuous(trans="log1p") +
  scale_y_continuous(trans="log1p")

ggplot(morph.df, aes(Holdfast, Weight_Holdfast, colour=source)) + 
  geom_point(alpha=0.5, shape=1) + 
  stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site)) +
  stat_smooth(method="lm", se=F, size=0.75, aes(group=source)) +
  scale_x_continuous("Holdfast diameter", trans="log1p") +
  scale_y_continuous(trans="log1p")

ggplot(morph.df, aes(Holdfast, Weight_Stipes, colour=source)) + 
  geom_point(alpha=0.5, shape=1) + 
  stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site)) +
  stat_smooth(method="lm", se=F, size=0.75, aes(group=source)) +
  scale_x_continuous("Holdfast diameter", trans="log1p") +
  scale_y_continuous(trans="log1p")

ggplot(morph.df, aes(Holdfast, Weight_Blades, colour=source)) + 
  geom_point(alpha=0.5, shape=1) + 
  stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site)) +
  stat_smooth(method="lm", se=F, size=0.75, aes(group=source)) +
  scale_x_continuous("Holdfast diameter", trans="log1p") +
  scale_y_continuous(trans="log1p")

# Larger holdfast = less blade-y by weight, more holdfast-y & stipe-y
ggplot(morph.df, aes(Holdfast, Weight_Holdfast/WeightTot, colour=source)) + 
  geom_point(alpha=0.5, shape=1) + 
  stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site)) +
  stat_smooth(method="lm", se=F, size=0.75, aes(group=source)) +
  scale_x_continuous("Holdfast diameter", trans="log1p") +
  scale_y_continuous("Holdfast weight / Total", limits=c(0,1))

ggplot(morph.df, aes(Holdfast, Weight_Stipes/WeightTot, colour=source)) + 
  geom_point(alpha=0.5, shape=1) + 
  stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site)) +
  stat_smooth(method="lm", se=F, size=0.75, aes(group=source)) +
  scale_x_continuous("Holdfast diameter", trans="log1p") +
  scale_y_continuous("Stipe weight / Total", limits=c(0,1))

ggplot(morph.df, aes(Holdfast, Weight_Blades/WeightTot, colour=source)) + 
  geom_point(alpha=0.5, shape=1) + 
  stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site)) +
  stat_smooth(method="lm", se=F, size=0.75, aes(group=source)) +
  scale_x_continuous("Holdfast diameter", trans="log1p") +
  scale_y_continuous("Blade weight / Total", limits=c(0,1))

ggplot(morph.df, aes(Holdfast, Height_Holdfast, colour=source)) + 
  geom_point(alpha=0.5, shape=1) + 
  stat_smooth(method="lm", se=F, size=0.25, aes(group=Site)) +
  scale_x_log10() + scale_y_log10()

ggplot(morph.df, aes(Holdfast, colour=factor(Frondosity))) + 
  geom_density() + scale_colour_viridis_d() +
  scale_x_log10() 

# Lots of NAs -- Cata's thesis should have this too
ggplot(morph.df, aes(Holdfast, colour=factor(Reproductive))) + 
  geom_density() + 
  scale_x_log10() 

2 Density

To come.

3 Growth & Survival dataset

Growth_Surv <- read_csv(glue("{data.dir}Growth&Survival.csv")) %>%
  mutate(Site=factor(Site),
         Management=factor(Management),
         Plant=factor(Plant),
         alive=as.numeric(Holdfast > 0),
         Date_Harvest=Date_harvest_dmy %>%
           str_replace("Dic", "01-12") %>%
           str_replace("Nov", "01-11") %>%
           str_replace("Ene", "01-01") %>%
           str_replace("Feb", "01-02") %>%
           dmy(),
         Date_Obs=Date_monitoring %>%
           str_replace("Abr", "01-04") %>%
           str_replace("Oct", "01-10") %>%
           str_replace("May", "01-05") %>%
           dmy()) %>%
  mutate(days_since_harvest=as.numeric(Date_Obs - Date_Harvest)) 

survey.i <- Growth_Surv %>%
  group_by(Site, Management, Patch, Monitoring) %>%
  summarise(Date_Harvest=min(Date_Harvest),
            Date_Obs=min(Date_Obs)) %>%
  arrange(Management, Site, Patch, Monitoring)
survey.i <- bind_rows(
  survey.i,
  survey.i %>% slice_head(n=1) %>%
    mutate(Monitoring=0,
           Date_Obs=Date_Harvest)
  ) %>%
  arrange(Management, Site, Patch, Monitoring) %>%
  mutate(days_since_harvest=as.numeric(Date_Obs - Date_Harvest),
         days_since_prevObs=as.numeric(Date_Obs - lag(Date_Obs)),
         days_since_prevObs=replace_na(days_since_prevObs, 0),
         days_to_nextObs=as.numeric(lead(Date_Obs) - Date_Obs))

Growth_Surv <- Growth_Surv %>%
  select(Site, Management, Patch, Monitoring, 
         Plant, Holdfast, Stipes, Total_Length) %>%
  left_join(., survey.i) %>%
  arrange(Site, Patch, Plant, Monitoring) %>%
  mutate(prevSurv=Monitoring-1) %>%
  arrange(Site, Patch, Plant, Monitoring) %>%
  group_by(Site, Patch, Plant) %>%
  mutate(doesSurvive=lead(Holdfast>0),
         nextHoldfast=lead(Holdfast),
         nextStipes=lead(Stipes),
         nextLength=lead(Total_Length),
         maxAge_days=cumsum(days_since_prevObs),
         maxAge_yrs=maxAge_days/365,
         surv_01=as.numeric(doesSurvive),
         logHoldfast=log(Holdfast)) %>%
  ungroup

4 Survival

We can link probability of survival to holdfast size x Management. With the magic of Hierarchical Bayesian modelling, we can annualize the estimated survival rates if we assume that the survival probability is constant throughout the year.

surv.df <- Growth_Surv %>% filter(!is.na(doesSurvive))

surv.df %>% group_by(Management) %>%
  summarise(nPlants=n_distinct(paste(Site, Patch, Plant)),
            nMortalities=sum(!doesSurvive),
            nSurvivals=sum(doesSurvive))

## # A tibble: 2 × 4
##   Management nPlants nMortalities nSurvivals
##   <fct>        <int>        <int>      <int>
## 1 OA              72           59         24
## 2 TURF            53           17         57

ggplot(surv.df, aes(Holdfast, as.numeric(doesSurvive))) +  
  geom_jitter(shape=1, height=0.025, width=0.25) +
  stat_smooth(method="glm", method.args=list(family="binomial")) +
  scale_x_continuous(trans="log", breaks=c(0.1, 1, 5, 20))

ggplot(surv.df, aes(Holdfast, as.numeric(doesSurvive), colour=Management)) +  
  geom_jitter(shape=1, height=0.025, width=0.25) +
  stat_smooth(method="glm", method.args=list(family="binomial")) +
  scale_x_continuous(trans="log", breaks=c(0.1, 1, 5, 20))

# This is asking a lot of the data, but we can get reasonable looking curves
# across management types. The outcome isn't all that sensitive to priors 
# at least within reason. Here's an option expecting greater variability among
# sites, then among plants, then among patches, and some nudged slopes.
library(brms)
out.s <- brm(bf(surv_01 ~ inv_logit(pSurvYr)^maxAge_yrs,
                pSurvYr ~ logHoldfast*Management + 
                  (1+logHoldfast*Management|Site/Patch/Plant),
                nl=T),
             data=surv.df, init=0, family="bernoulli", chains=4, cores=4,
             prior=c(prior(normal(0, 1), nlpar="pSurvYr"),
                     prior(normal(-1, 1), nlpar="pSurvYr", coef="Intercept"),
                     prior(normal(1, 1), nlpar="pSurvYr", coef="ManagementTURF"),
                     prior(normal(0, 0.5), nlpar="pSurvYr", class="sd", lb=0, 
                           group="Site"),
                     prior(normal(0, 0.1), nlpar="pSurvYr", class="sd", lb=0,
                           group="Site:Patch"),
                     prior(normal(0, 0.25), nlpar="pSurvYr", class="sd", lb=0,
                           group="Site:Patch:Plant")))

fitted.df <- expand_grid(logHoldfast=seq(-2.5,3.5,length.out=100),
            Management=levels(surv.df$Management)) %>%
  mutate(maxAge_yrs=1,
         Holdfast=exp(logHoldfast),
         Site=NA, Patch=NA, Plant=NA) %>%
  bind_cols(posterior_epred(out.s, newdata=., nlpar="pSurvYr", 
                            ndraws=4000, re_formula=NA) %>%
                        t %>% as_tibble) %>%
  pivot_longer(starts_with("V"), names_to="iter", values_to="logit_p") %>%
  group_by(logHoldfast, Holdfast, Management) %>%
  summarise(logit_p_mn=mean(logit_p),
            logit_p_lo=tidybayes::qi(logit_p, 0.9)[,1],
            logit_p_hi=tidybayes::qi(logit_p, 0.9)[,2]) %>%
  mutate(across(starts_with("logit_p"), inv_logit_scaled, 
                .names="{str_sub({.col}, 7, -1)}"))
fitted.df %>%
  ggplot(aes(Holdfast, p_mn, colour=Management, fill=Management)) + 
  geom_ribbon(aes(ymin=p_lo, ymax=p_hi), alpha=0.5, colour=NA) +
  geom_line() +
  geom_jitter(data=surv.df, aes(y=surv_01),
              shape=1, height=0.025, width=0.05) +
  scale_x_continuous(trans="log", breaks=c(0.1, 1, 5, 20)) +
  labs(x="Holdfast diameter (cm)", y="Annual survival probability")

fitted.df %>%
  ggplot(aes(logHoldfast, logit_p_mn, colour=Management, fill=Management)) + 
  geom_ribbon(aes(ymin=logit_p_lo, ymax=logit_p_hi), alpha=0.5, colour=NA) +
  geom_line() +
  geom_jitter(data=surv.df, aes(y=logit_scaled(pmax(pmin(surv_01, 0.999), 0.001))),
              shape=1, height=0.025, width=0.05) +
  labs(x="log holdfast diameter (cm)", y="Annual survival probability")

5 Growth

We can link the holdfast size at t+1 to holdfast size at t x Management. Like with mortality, we can assume a constant daily growth rate to account for differences in the number of days between surveys and calculate annualized growth rates by holdfast diameter.

grow.df <- Growth_Surv %>% 
  filter(doesSurvive) %>%
  mutate(annHoldfastGrowth=(nextHoldfast-Holdfast)/days_to_nextObs*365,
         annStipeGrowth=(nextStipes-Stipes)/days_to_nextObs*365,
         annLengthGrowth=(nextLength-Total_Length)/days_to_nextObs*365,
         nextHoldfastAnnual=Holdfast+annHoldfastGrowth,
         logNextHoldfastAnnual=log(nextHoldfastAnnual))

ggplot(grow.df, aes(Holdfast, Holdfast+annHoldfastGrowth)) + 
  geom_abline(colour="grey") +
  geom_point() + stat_smooth(method="lm", formula=y~x+I(x^2)) +
  scale_x_log10() + scale_y_log10()

ggplot(grow.df, aes(Holdfast, Holdfast+annHoldfastGrowth, colour=Management)) +
  geom_abline(colour="grey") +
  geom_point() + 
  stat_smooth(aes(group=Site), method="lm", formula=y~x+I(x^2), se=F, size=0.25) +
  stat_smooth(method="lm", formula=y~x+I(x^2)) +
  scale_x_log10() + scale_y_log10()

ggplot(grow.df, aes(Stipes, Stipes+annStipeGrowth)) +
  geom_abline(colour="grey") +
  geom_point() + stat_smooth(method="lm", formula=y~x+I(x^2)) 

ggplot(grow.df, aes(Stipes, Stipes+annStipeGrowth, colour=Management)) + 
  geom_abline(colour="grey") +
  geom_point() + 
  stat_smooth(aes(group=Site), method="lm", formula=y~x+I(x^2), se=F, size=0.25) +
  stat_smooth(method="lm", formula=y~x+I(x^2)) 

ggplot(grow.df, aes(Total_Length, Total_Length+annLengthGrowth)) + 
  geom_abline(colour="grey") +
  geom_point() + stat_smooth(method="lm", formula=y~x+I(x^2))

ggplot(grow.df, aes(Total_Length, Total_Length+annLengthGrowth, colour=Management)) +
  geom_abline(colour="grey") +
  geom_point() + 
  stat_smooth(aes(group=Site), method="lm", formula=y~x+I(x^2), se=F, size=0.25) +
  stat_smooth(method="lm", formula=y~x+I(x^2)) 

# log-log scale, heteroscedastic
out.g <- brm(bf(logNextHoldfastAnnual ~ logHoldfast*Management + 
                  (1+logHoldfast|Site/Patch/Plant),
                sigma ~ logHoldfast*Management),
             family="gaussian", data=grow.df, cores=4,
             prior=c(prior(normal(0, 0.1), class="sd", lb=0, group="Site"),
                     prior(normal(0, 0.1), class="sd", lb=0, group="Site:Patch"),
                     prior(normal(0, 0.1), class="sd", lb=0, group="Site:Patch:Plant"),
                     prior(normal(0, 1), class="b"),
                     prior(normal(0, 1), class="b", dpar="sigma"),
                     prior(normal(2, 1), class="Intercept"),
                     prior(normal(0, 1), class="Intercept", dpar="sigma")))

fitted.df <- expand_grid(Holdfast=seq(0.1,44,length.out=100),
                         Management=levels(surv.df$Management)) %>%
  mutate(yrs_to_next=1,
         logHoldfast=log(Holdfast),
         Site=NA, Patch=NA, Plant=NA)
preds <- posterior_epred(out.g, newdata=fitted.df, re.form=NA)
# preds <- posterior_predict(out.g, newdata=fitted.df, re.form=NA)
fitted.df <- fitted.df %>%
  mutate(pred_mn=colMeans(preds),
         pred_lo=apply(preds, 2, function(x) quantile(x, probs=0.1)),
         pred_hi=apply(preds, 2, function(x) quantile(x, probs=0.9)),
         pred_emn=colMeans(exp(preds)),
         pred_elo=apply(exp(preds), 2, function(x) quantile(x, probs=0.1)),
         pred_ehi=apply(exp(preds), 2, function(x) quantile(x, probs=0.9)))
fitted.df %>%
  ggplot(aes(logHoldfast, pred_mn, colour=Management, fill=Management)) + 
  geom_abline() +
  geom_ribbon(aes(ymin=pred_lo, ymax=pred_hi), alpha=0.5, colour=NA) +
  geom_line() +
  geom_point(data=grow.df, aes(y=log(Holdfast+annHoldfastGrowth)), shape=1) +
  labs(x="Holdfast diameter (t)", y="Holdfast diameter (t+1)")

fitted.df %>%
  ggplot(aes(Holdfast, pred_emn, colour=Management, fill=Management)) + 
  geom_abline() +
  geom_ribbon(aes(ymin=pred_elo, ymax=pred_ehi), alpha=0.5, colour=NA) +
  geom_line() +
  geom_point(data=grow.df, aes(y=Holdfast+annHoldfastGrowth), shape=1) +
  labs(x="Holdfast diameter (t)", y="Holdfast diameter (t+1)") +
  xlim(0, 50) + ylim(0, 50)

6 Recruitment

We can estimate size distributions of recruits by management, and estimate the number of recruits as a function of the number or size of non-recruits in the plot. TODO: Account for the differences in time between monitoring periods to annualize somehow.

rcr.df <- Growth_Surv %>% 
  select(Site, Management, Patch, Monitoring, 
         Plant, Holdfast, Stipes, Total_Length) %>%
  left_join(., survey.i) %>%
  arrange(Site, Patch, Plant, Monitoring) %>%
  mutate(prevSurv=Monitoring-1) %>%
  arrange(Site, Patch, Plant, Monitoring) %>%
  group_by(Site, Patch, Plant) %>%
  mutate(newRecruit=is.na(lag(Holdfast)),
         maxAge_days=cumsum(days_since_prevObs),
         maxAge_yrs=maxAge_days/365) %>%
  group_by(Site, Patch, Monitoring) %>%
  mutate(nPlants=sum(Holdfast>0),
         totHoldfast=sum(Holdfast)) %>%
  ungroup %>%
  filter(newRecruit==TRUE) %>%
  filter(maxAge_yrs < 1.5) 


# Some bias in measurements: preference for 0.2, 10 cm
table(rcr.df$Holdfast)

## 
## 0.1 0.2 0.6 0.8   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16 
##   8  60   1   1   8   6   6   5   3   6   1   3   1  10   1   6   1   1   3   3 
##  18  20  22  23  24  28  30 
##   2   1   2   1   1   1   2

ggplot(rcr.df, aes(Holdfast)) + geom_density() + 
  ggtitle("Recruit size distr.") +
  scale_x_continuous(trans="log1p")

ggplot(rcr.df, aes(Holdfast, colour=Management)) + 
  geom_density(aes(group=Site), size=0.25) + 
  geom_density(size=1) +
  ggtitle("Recruit size distr.") +
  scale_x_continuous(trans="log1p")

ggplot(rcr.df, aes(maxAge_yrs)) + geom_histogram() +
  ggtitle("Recruit max possible age")

ggplot(rcr.df, aes(Holdfast)) + geom_histogram() +
  scale_x_continuous(trans="log1p") + facet_grid(Management~Monitoring)

ggplot(rcr.df, aes(factor(Monitoring), fill=Management)) + 
  geom_bar(position="dodge") + 
  ggtitle("Number of recruits decreases with time")

ggplot(rcr.df, aes(Management, fill=factor(Monitoring))) + 
  geom_bar(position="fill") + 
  ggtitle("Number of recruits decreases with time, maybe more strongly in TURF")

rcr.df %>% group_by(Site, Management, Patch, Monitoring, nPlants) %>%
  summarise(nRecruits=sum(newRecruit)) %>%
  ggplot(aes(nPlants-nRecruits, nRecruits)) + 
  geom_jitter(width=0.2, height=0.2) +
  stat_smooth(method="lm", formula=y~log1p(x)) +
  scale_y_continuous(trans="log1p") +
  xlab("Number of non-recruits")

rcr.df %>% group_by(Site, Management, Patch, Monitoring, nPlants) %>%
  summarise(nRecruits=sum(newRecruit)) %>%
  ggplot(aes(nPlants-nRecruits, nRecruits, colour=Management)) + 
  geom_jitter(width=0.2, height=0.2) +
  stat_smooth(aes(group=Site), method="lm", formula=y~log1p(x), se=F, size=0.25) +
  stat_smooth(method="lm", formula=y~log1p(x), size=1) +
  scale_y_continuous(trans="log1p") +
  xlab("Number of non-recruits")

rcr.df %>% group_by(Site, Management, Patch, Monitoring, totHoldfast) %>%
  summarise(nRecruits=sum(newRecruit),
            HoldfastRcr=sum(newRecruit*Holdfast)) %>%
  ggplot(aes(totHoldfast-HoldfastRcr, nRecruits)) + 
  geom_jitter(width=0.2, height=0.2) +
  stat_smooth(method="lm", formula=y~log1p(x)) +
  scale_y_continuous(trans="log1p") +
  xlab("Total holdfast diameter of non-recruits")

rcr.df %>% group_by(Site, Management, Patch, Monitoring, totHoldfast) %>%
  summarise(nRecruits=sum(newRecruit),
            HoldfastRcr=sum(newRecruit*Holdfast)) %>%
  ggplot(aes(totHoldfast-HoldfastRcr, nRecruits, colour=Management)) + geom_jitter() +
  stat_smooth(aes(group=Site), method="lm", formula=y~log1p(x), se=F, size=0.25) +
  stat_smooth(method="lm", formula=y~log1p(x), size=1) +
  scale_y_continuous(trans="log1p") +
  xlab("Total holdfast diameter of non-recruits")

# The recruit size distribution seems to change with time-since-harvest.
# Mainly small recruits at the start, then larger recruits at TURF sites as 
# populations become established... These should only be *new* recruits since
# the last survey. Differences in thinning / establishment by management? OA
# tend to have more and smaller recruits after the initial pulse. 

# How does this interact with higher survival rates and faster growth rates in TURF?
# Lower establishment, but better survival and growth once established?
ggplot(rcr.df, aes(days_since_harvest, Holdfast)) + 
  geom_jitter(width=10, height=0.1) + 
  stat_smooth(span=1.5) +
  ggtitle("Recruit size vs. time since clearance") +
  scale_y_continuous(trans="log1p")

ggplot(rcr.df, aes(days_since_harvest, Holdfast, colour=Management)) + 
  geom_jitter(width=10, height=0.1) + 
  stat_smooth(aes(group=Site), span=1.5, se=F, size=0.25) +
  stat_smooth(span=1.5) +
  ggtitle("Recruit size vs. time since clearance") +
  scale_y_continuous(trans="log1p")

ggplot(rcr.df, aes(Holdfast, colour=factor(Monitoring))) + 
  geom_density(aes(group=paste(Site, Monitoring)), adjust=1.5) + 
  ggtitle("Recruit size vs. time since clearance") +
  scale_y_continuous(trans="log1p") +
  facet_wrap(~Management)

rcr.df %>% group_by(Site, Management, Patch, Monitoring, days_since_harvest) %>%
  summarise(nRecruits=sum(newRecruit)) %>%
  ggplot(aes(days_since_harvest, nRecruits)) + 
  geom_jitter(width=10, height=0.1) + 
  stat_smooth(span=1.5) +
  ggtitle("Number of recruits vs. time since clearance") +
  scale_y_continuous(trans="log1p")

rcr.df %>% group_by(Site, Management, Patch, Monitoring, days_since_harvest) %>%
  summarise(nRecruits=sum(newRecruit)) %>%
  ggplot(aes(days_since_harvest, nRecruits, colour=Management)) + 
  geom_jitter(width=10, height=0.1) + 
  stat_smooth(aes(group=Site), span=1.5, se=F, size=0.25) +
  stat_smooth(span=1.5) +
  ggtitle("Number of recruits vs. time since clearance") +
  scale_y_continuous(trans="log1p")

7 Reproduction

We can fit a nice curve predicting the probability of reproductive blades given holdfast diameter, but do we have any data linking reproductive blades to recruits? That is, what we need is \(F(z',z) \sim p_{repro}(z) f_{spores}(z) p_{estab} f_{rcrSize}(z')\). We have an estimate for \(p_{repro}(z)\), which is the probability of having reproductive blades given a holdfast size \(z\). We have part of \(f_{rcrSize}(z')\), which is the expected recruit size distribution at time \(t+1\). What we’re missing is how many spores might be produced by each reproductive plant, and how likely those spores are to establish.

We could potentially lump several processes and calculate directly the number of recruits from the data above based on a presumed maximum established recruits, with the actual as a function of \(N\).

morph.df %>%
  filter(!is.na(Reproductive)) %>% 
  droplevels() %>%
  ggplot(aes(Holdfast, Reproductive)) + 
  geom_jitter(shape=1, height=0.025, width=0.5) +
  stat_smooth(method="glm", method.args=list(family="binomial")) +
  scale_x_continuous(trans="log1p")

morph.df %>%
  filter(!is.na(Reproductive)) %>% 
  droplevels() %>%
  ggplot(aes(Holdfast, Reproductive, colour=Management)) + 
  geom_jitter(shape=1, height=0.025, width=0.5) +
  stat_smooth(method="glm", method.args=list(family="binomial")) +
  scale_x_continuous(trans="log1p")